CL_M <-
function(x,y,s,mu,A,c,inds){
	n <- length(inds)
	k <- length(x)
	out <- matrix(0,k,k)
	for(i in 1:n){
		ki <- length(inds[[i]])
		w1bar <- x[inds[[i]]] + mu[i]*s[inds[[i]]]
		w2bar <- mu[i]*x[inds[[i]]] + s[inds[[i]]]
		wbar <- jordan_sqrt(jordan_square(w1bar) + jordan_square(w2bar) + 2*(mu[i]^2)*jordan_identity(ki))
		Lwinv <- solve(jordan_sym(wbar))
		
		out[inds[[i]],inds[[i]]] <- (1+mu[i])*diag(ki) - Lwinv %*% (jordan_sym(w1bar) + mu[i]*jordan_sym(w2bar))	
	}
	return(out)
}

